Aplicações de Modelos Lineares Generalizados

Pacotes

library(tidyverse)
library(psych)
library(plotly)
library(GGally)
library(car)
library(gridExtra)
library(cowplot)
library(DescTools)

I. Modelo Linear Normal

Para a aplicação a seguir, consideraremos o modelo de regressão normal linear

\[ y_i = \beta_1 + \beta_2x_{2i} + ... + \beta_px_{pi} + \epsilon_i, \quad i = 1,...,n \]

em que os erros \(e_i\) são variáveis aleatórias independentes, normalmente distribuídas, de média zero e variância \(\sigma^2\) constante.

1. Preço de Venda de Imóveis

Dados

Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 23, página 112). Dados disponíveis em: <https://www.ime.usp.br/~giapaula/textoregressao.htm>, arquivo imoveis.dat.

Neste exemplo, vamos modelar o preço de venda de imóveis a partir de dados relativos a uma amostra de 27 imóveis. As variáveis do conjunto de dados são:

• imposto: imposto do imóvel (em US$ 100);

• areaT: área do terreno (em 1.000 pés quadrados);

• areaC: área construída (em 1.000 pés quadrados);

• idade: idade da residência (em anos);

• preco: preço de venda do imóvel (em US$ 1.000).

Sendo assim, o nosso objetivo é encontrar o melhor ajuste linear que explique e quantifique a variável preço de venda a partir das demais.

# Obtendo os dados
dados = read.table("dados/imoveis.dat")
colnames(dados) = c("imposto", "areaT", "areaC", "idade", "preco")
attach(dados)

# Algumas observações dos dados
head(dados)

##   imposto areaT areaC idade preco
## 1  4.9176 3.472 0.998    42  25.9
## 2  5.0208 3.531 1.500    62  29.5
## 3  4.5429 2.275 1.175    40  27.9
## 4  4.5573 4.050 1.232    54  25.9
## 5  5.0597 4.455 1.121    42  29.9
## 6  3.8910 4.455 0.988    56  29.9

# Medidas descritivas
describe(dados)

##         vars  n  mean    sd median trimmed   mad   min   max range  skew
## imposto    1 27  7.24  2.88   6.09    6.84  2.15  3.89 15.42 11.53  1.40
## areaT      2 27  6.35  2.40   5.85    6.22  2.07  2.28 12.80 10.53  0.68
## areaC      3 27  1.51  0.56   1.49    1.41  0.39  0.98  3.42  2.44  2.02
## idade      4 27 36.48 14.05  40.00   36.96 14.83  3.00 62.00 59.00 -0.34
## preco      5 27 38.50 14.31  36.90   35.65  8.90 25.90 84.90 59.00  2.25
##         kurtosis   se
## imposto     1.34 0.55
## areaT       0.01 0.46
## areaC       4.08 0.11
## idade      -0.54 2.70
## preco       4.63 2.75

Análise Descritiva

Box-plots

plot_ly(type = "box") %>%
  add_trace(y = imposto, name = "imposto") %>% 
  add_trace(y = areaT, name = "área do terreno") %>% 
  add_trace(y = areaC, name = "área construída") %>%
  add_trace(y = idade, name = "idade do imóvel") %>%
  add_trace(y = preco, name = "preço de venda")

Em um breve resumo descritivo dos dados, observamos que as altas variâncias das variáveis idade e preço destacam-se das demais. Com uma idade mediana de 40 anos, temos um imóvel de apenas 3 anos de idade e, com um preço mediano de US$ 36.900, temos dois imóveis bem mais caros nos valores de US$ 82.900 e US$ 84.900.

Os gráficos de box-plot nos mostram que os dados possuem outliers e que as variáveis possuem distribuições assimétricas. Focando na variável alvo, preço de venda, os pontos destoantes são justamente os dois imóveis mais caros.

Densidades, Dispersões e Correlações

g = ggpairs(dados, aes(color = I("slategray"), fill = I("slategray")),
            lower = list(continuous = wrap("smooth",col="black")),
            diag=list(continuous = wrap("densityDiag",alpha=0.5,size=1)))
ggplotly(g)

No gráfico acima, as curvas representadas na diagonal principal deixam mais evidente a constatação de assimetria nas distribuições observada nos box-plots. A variável preço possui correlações positivas fortes com as variáveis imposto, área do terreno e área construída, mas uma correlação negativa fraca com a variável idade. Em outras palavras, quanto maiores os valores de imposto, área construída e área do terreno, maior o preço de venda do imóvel.

Podemos observar também que há correlações relevantes entre as variáveis explicativas imposto, área construída e área do terreno. Isso nos dá indícios de multicolinearidade e, assim, uma possível redundância de informação passada por essas variáveis. Então, talvez um modelo completo com todas as variáveis não seja o mais adequado.

Modelo Completo

Nesse primeiro ajuste, consideraremos o modelo completo com todas as variáveis disponíveis.

\[ preço_i = \beta_0 + \beta_1imposto_i + \beta_2areaC_i + \beta_3areaT_i + \beta_4idade_i \]

model1 = lm(preco~.,dados)
summary(model1)

## 
## Call:
## lm(formula = preco ~ ., data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9685 -2.7152  0.2663  2.1374  6.9872 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.43587    4.09233   0.595  0.55776    
## imposto      2.07823    0.55296   3.758  0.00109 ** 
## areaT        0.23238    0.50658   0.459  0.65094    
## areaC       13.97381    2.90663   4.808  8.4e-05 ***
## idade       -0.04376    0.06628  -0.660  0.51594    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.077 on 22 degrees of freedom
## Multiple R-squared:  0.9313, Adjusted R-squared:  0.9188 
## F-statistic: 74.56 on 4 and 22 DF,  p-value: 1.809e-12

Verificamos que, para o modelo completo, as variáveis imposto e área construída são as únicas significativas e obtivemos um \(R^2\) ajustado de 0.92, muito alto, indicando um possível bom ajuste aos dados amostrais. Porém, o \(R^2\) por si só não nos garante que o modelo seja o mais adequado.

Seleção de variáveis

Aplicaremos o critério de informação de Akaike (AIC) ao modelo completo, para selecionar as variáveis que deverão permanecer.

MASS::stepAIC(model1)

## Start:  AIC=80.36
## preco ~ imposto + areaT + areaC + idade
## 
##           Df Sum of Sq    RSS    AIC
## - areaT    1      3.50 369.14 78.614
## - idade    1      7.25 372.88 78.887
## <none>                 365.64 80.357
## - imposto  1    234.76 600.40 91.748
## - areaC    1    384.13 749.77 97.746
## 
## Step:  AIC=78.61
## preco ~ imposto + areaC + idade
## 
##           Df Sum of Sq    RSS    AIC
## - idade    1     11.51 380.64 77.443
## <none>                 369.14 78.614
## - imposto  1    246.10 615.24 90.407
## - areaC    1    489.34 858.47 99.402
## 
## Step:  AIC=77.44
## preco ~ imposto + areaC
## 
##           Df Sum of Sq    RSS    AIC
## <none>                 380.64 77.443
## - imposto  1    350.01 730.65 93.049
## - areaC    1    483.16 863.80 97.569

## 
## Call:
## lm(formula = preco ~ imposto + areaC, data = dados)
## 
## Coefficients:
## (Intercept)      imposto        areaC  
##      0.7903       2.2971      13.9333

O AIC nos retornou as variáveis imposto e área construída como aquelas que deverão ser mantidas no modelo, o qual chamaremos de modelo parcimonioso.

Modelo Parcimonioso

model2 = lm(preco~imposto+areaC,dados)
summary(model2)

## 
## Call:
## lm(formula = preco ~ imposto + areaC, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0177 -3.3169  0.0958  2.4382  7.0949 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7903     2.2794   0.347    0.732    
## imposto       2.2971     0.4890   4.698 8.96e-05 ***
## areaC        13.9333     2.5244   5.519 1.12e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.982 on 24 degrees of freedom
## Multiple R-squared:  0.9285, Adjusted R-squared:  0.9225 
## F-statistic: 155.8 on 2 and 24 DF,  p-value: 1.791e-14

E, assim, obtemos um modelo em que as variáveis imposto e área construída são altamente significativas.

Porém, na análise descritiva, observamos que as variáveis área construída e área do terreno possuem uma correlação considerável. Então, podemos ainda testar um segundo modelo parcimonioso considerando a variável área do terreno ao invés da área construída.

model3 = lm(preco~imposto+areaT,dados)
summary(model3)

## 
## Call:
## lm(formula = preco ~ imposto + areaT, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.621  -3.070   0.446   3.167  10.610 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.1614     3.3019   0.957   0.3479    
## imposto       3.9126     0.5293   7.391 1.24e-07 ***
## areaT         1.1017     0.6347   1.736   0.0954 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.655 on 24 degrees of freedom
## Multiple R-squared:  0.8558, Adjusted R-squared:  0.8438 
## F-statistic: 71.22 on 2 and 24 DF,  p-value: 8.081e-11

E, assim, obtemos um modelo em que a variável área do terreno passa de não significativa para pouco significativa ao nível de 10% de significância. Já a variável imposto continua bastante significativa ao nível de 1% de significância.

Diagnóstico

Agora, faremos a análise de diagnóstico dos dois modelos parcimoniosos propostos, para verificar suas adequações.

Envelope

par(mfrow=c(1,2))
qqPlot(model2, pch = 16, main = "preco ~ imposto + areaC")

## [1]  6 10

qqPlot(model3, pch = 16, main = "preco ~ imposto + areaT")

## [1]  9 27

Teste de Normalidade

\(H_0\): Os resíduos têm distribuição normal.

# Shapiro-Wilk
shapiro.test(model2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model2$residuals
## W = 0.96558, p-value = 0.4903

shapiro.test(model3$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model3$residuals
## W = 0.9793, p-value = 0.8462

Ao nível de 5% de significância, não rejeitamos a hipótese de que os resíduos tenham distribuição normal. Nos gráficos de envelope para o ajuste normal, os pontos estão dentro das bandas, exceto pela observação 27 no modelo com área do terreno. Porém, a presença de leves ondulações pode ser indício de variância não constante. Vamos verificar se de fato isso ocorre no gráfico de resíduos por valores ajustados.

Resíduos x Valores ajustados

residualPlot(model2, type = "rstudent", id = TRUE,
             ylab = "resíduo studentizado", xlab = "valor ajustado",
             main = "preco ~ imposto + areaC", pch = 16, quadratic = FALSE)

residualPlot(model3, type = "rstudent", id = TRUE,
             ylab = "resíduo studentizado", xlab = "valor ajustado",
             main = "preco ~ imposto + areaT", pch = 16, quadratic = FALSE)

Embora não haja violação do pressuposto de normalidade dos resíduos, o gráfico de resíduos por valores ajustados nos dá indícios de heterocedasticidade, isto é, variância não constante, como havíamos suposto a partir do gráfico de envelope. Porém, esse comportamento pode estar sendo influenciado pelas observações destoantes vistas na etapa descritiva. Além disso, podemos perceber no gráfico de resíduos por valores ajustados a presença de pontos aberrantes, isto é, muito além do intervalo [-2,2]. No caso do primeiro modelo parcimonioso, o ponto aberrante é a observação 10, enquanto que no segundo modelo parcimonioso são as observações 9 e 27. Sendo assim, a seguir, vamos verificar se de fato esses pontos estão influenciando na qualidade dos ajustes.

Pontos Influentes

Partiremos agora para a investigação de possíveis pontos influentes, que podem estar atrapalhando a qualidade do ajuste.

Medidas de influência:

• DFBetas (dfb): estatísticas que indicam o efeito da remoção de cada observação sobre as estimativas dos parâmetros do modelo;

• DFFit (dffit) e Cook’s D (cook.d): são estatísticas que indicam o efeito da remoção de cada observação sobre os valores ajustados do modelo;

• COVRATIO (cov.r): estatística que indica o efeito da remoção de cada observação sobre a matriz de covariâncias do modelo, em outras palavras, mede a alteração na precisão das estimativas dos parâmetros do modelo;

• HAT (hat): diagonal da matriz de projeção (\(H = X(X'X)^{-1}X'\)) da solução dos mínimos quadrados, é a métrica de alavancagem.

inf = influence.measures(model2)
summary(inf)

## Potentially influential observations of
##   lm(formula = preco ~ imposto + areaC, data = dados) :
## 
##    dfb.1_  dfb.imps dfb.areC dffit   cov.r   cook.d  hat    
## 9  -0.28    0.00     0.19     0.35    2.18_*  0.04    0.49_*
## 10 -1.25_*  0.29     0.59     1.62_*  0.87    0.73    0.32  
## 27 -0.13   -2.04_*   1.86_*  -2.12_*  2.04_*  1.39_*  0.61_*

inf = influence.measures(model3)
summary(inf)

## Potentially influential observations of
##   lm(formula = preco ~ imposto + areaT, data = dados) :
## 
##    dfb.1_  dfb.imps dfb.areT dffit   cov.r   cook.d  hat    
## 9  -1.15_*  1.65_*  -0.46     2.00_*  0.81    1.07_*  0.37_*
## 10 -1.21_*  0.47     0.68     1.54_*  1.02    0.69    0.35_*
## 27  0.09   -2.83_*   2.36_*  -3.05_*  0.33_*  1.85_*  0.35_*

Baseando-se nas medidas de influência acima, que verificam o efeito da remoção de cada uma das observações individualmente, temos que os possíveis pontos de influência são as observações 9 e 27, para ambos os modelos parcimoniosos, e a observação 10, para o segundo. Percebemos também que a observação 27 é a mais crítica, por impactar nas estimativas dos dois parâmetros, nos valores ajustados e na variância do modelo. Também podemos constatar isso graficamente, como se segue a baixo.

# DISTANCIA DE COOK

plot(model2, which=4, lwd=5, main="Distância de Cook x Observação", caption=F,
     sub.caption = "preco ~ imposto + areaC")

plot(model3, which=4, lwd=5, main="Distância de Cook x Observação", caption=F,
     sub.caption = "preco ~ imposto + areaT")

# ALAVANCAGEM

cut = 2*3/27 # 2*p/n
plot(hatvalues(model2), pch = 16,
     xlab = "observação", ylab = "medida h", main = "Alavancagem",
     sub = "preco ~ imposto + areaC")
abline(h = cut, lty = 2, lwd = 2)
text(hatvalues(model2) * as.integer(hatvalues(model2) > cut),
     cex=0.8, p=1, offset=0.3)

cut = 2*3/27 # 2*p/n
plot(hatvalues(model3), pch = 16,
     xlab = "observação", ylab = "medida h", main = "Alavancagem",
     sub = "preco ~ imposto + areaT")
abline(h = cut, lty = 2, lwd = 2)
text(hatvalues(model3) * as.integer(hatvalues(model3) > cut),
     cex=0.8, p=1, offset=0.3)

As observações 9, 10 e 27 são referentes aos dados:

dados[c(9,10,27),c("imposto","areaC","areaT","preco")]

##    imposto areaC areaT preco
## 9  15.4202  3.42   9.8  84.9
## 10 14.4598  3.00  12.8  82.9
## 27 12.0000  1.20   5.0  41.0

Os imóveis 9 e 10 possuem os valores de imposto mais elevados e as maiores áreas de terreno e de construção. De acordo com a relação linear identificada na etapa descritiva e a significância dessas variáveis para o modelo preditivo, podemos dizer que essas características justificam um alto preço de venda. Porém, os valores atribuídos aos imóveis 9 e 10 são muito mais elevados que os esperados pelo modelo proposto a imóveis com os mesmos atributos. Já a observação 27 causa estranheza por ter um valor de imposto tão alto e próximo aos dos imóveis 9 e 10, mesmo tendo menos da metade de suas áreas de terreno e de construção e custando pouco menos da metade que o imóvel 10.

Impacto das remoções das observações

# Compara os coeficientes e seus respectivos p-valores de um modelo
# após a remoção de uma observação
impacto = function(m1, m2){
  imp = (coef(m1) - coef(m2))/coef(m1) * 100
  impacto = data.frame(coef_com_i = coef(m1),
                       coef_sem_i = coef(m2),
                       impacto = imp,
                       pval_com_i = summary(m1)$coefficients[,4],
                       pval_sem_i = summary(m2)$coefficients[,4])
  return(impacto)
}

• Observação i=9:

  model2_sem9 = lm(preco~imposto+areaC,dados,subset=-9)
  impacto(model2,model2_sem9)

  ##             coef_com_i coef_sem_i      impacto   pval_com_i   pval_sem_i
  ## (Intercept)  0.7903027   1.433603 -81.39918902 7.318305e-01 0.6304908219
  ## imposto      2.2970580   2.297773  -0.03113799 8.955810e-05 0.0001221436
  ## areaC       13.9332510  13.454944   3.43284823 1.122859e-05 0.0001144541

  model3_sem9 = lm(preco~imposto+areaT,dados,subset=-9)
  impacto(model3,model3_sem9)

  ##             coef_com_i coef_sem_i    impacto   pval_com_i   pval_sem_i
  ## (Intercept)   3.161412   6.558508 -107.45502 3.478823e-01 5.383568e-02
  ## imposto       3.912561   3.131158   19.97165 1.243498e-07 1.077092e-05
  ## areaT         1.101699   1.360792  -23.51762 9.543188e-02 2.722222e-02

  A remoção da observação 9 não teve impactos relevantes nos coeficientes dos dois modelos parcimoniosos e também não causou mudança inferencial, uma vez que as variáveis imposto, área construída e área do terreno permaneceram igualmente significativas nos respectivos modelos;

• Observação i=10:

  model2_sem10 = lm(preco~imposto+areaC,dados,subset=-10)
  impacto(model2,model2_sem10)

  ##             coef_com_i coef_sem_i     impacto   pval_com_i   pval_sem_i
  ## (Intercept)  0.7903027   3.398847 -330.069021 7.318305e-01 1.640664e-01
  ## imposto      2.2970580   2.167849    5.624971 8.955810e-05 7.657851e-05
  ## areaC       13.9332510  12.571459    9.773685 1.122859e-05 2.390378e-05

  model3_sem10 = lm(preco~imposto+areaT,dados,subset=-10)
  impacto(model3,model3_sem10)

  ##             coef_com_i coef_sem_i     impacto   pval_com_i   pval_sem_i
  ## (Intercept)   3.161412  6.8915170 -117.988556 3.478823e-01 6.463495e-02
  ## imposto       3.912561  3.6816714    5.901246 1.243498e-07 2.123188e-07
  ## areaT         1.101699  0.6967438   36.757341 9.543188e-02 2.749152e-01

  A remoção da observação 10 não causou impactos relevantes nos coeficientes do primeiro modelo parcimonioso e também não ocasionou mudança inferencial, tendo as variáveis imposto e área construída mantido seus níveis de significância. Já no segundo modelo parcimonioso, a remoção da observação 10 ocasionou uma mudança inferencial, uma vez que a variável área de terreno perde sua significância;

• Observação i=27:

  model2_sem27 = lm(preco~imposto+areaC,dados,subset=-27)
  impacto(model2,model2_sem27)

  ##             coef_com_i coef_sem_i   impacto   pval_com_i   pval_sem_i
  ## (Intercept)  0.7903027   1.068075 -35.14754 7.318305e-01 0.6320276693
  ## imposto      2.2970580   3.255657 -41.73159 8.955810e-05 0.0001915344
  ## areaC       13.9332510   9.412139  32.44836 1.122859e-05 0.0155595300

  model3_sem27 = lm(preco~imposto+areaT,dados,subset=-27)
  impacto(model3,model3_sem27)

  ##             coef_com_i coef_sem_i    impacto   pval_com_i   pval_sem_i
  ## (Intercept)   3.161412  2.9423870   6.928084 3.478823e-01 2.603227e-01
  ## imposto       3.912561  5.0694489 -29.568550 1.243498e-07 4.765596e-10
  ## areaT         1.101699 -0.0528502 104.797154 9.543188e-02 9.260592e-01

  A remoção da observação 27 impactou em mudança inferencial nos dois modelos parcimoniosos. Enquanto que no primeiro, a significância da variável área construída é reduzida, no segundo, a variável área do terreno perde completamente sua significância.

Como o segundo modelo parcimonioso não nos trouxe nenhuma melhoria para o ajuste e as conclusões sobre os pontos de influência não favoreceram a permanência da variável área do terreno, optaremos pelo primeiro modelo parcimonioso com as variáveis imposto e área construída sugerido inicialmente pelo critério de seleção de Akaike.

A partir dos gráficos de resíduos por valores ajustados, de distância de Cook e de alavancagem, podemos fazer as seguintes classificações:

• Ponto aberrante: observação 10;

• Pontos de alavanca: observações 9, 10 e 27;

• Ponto influente: observação 27.

Interpretação

Então, concluímos que

\[ preço_i = 0.79 + 2.30*imposto_i + 13.93*areaC_i \quad , \]

onde:

• A cada US$ 100,00 acrescidos no valor de imposto, aumenta-se, em média, US$ 2.300,00 no preço de venda;

• A cada 1.000 metros quadrados de área construída acrescidos, aumenta-se, em média, US$ 13.930,00 no preço de venda.

II. Modelo Gama

1. Tempo de Vida de Pacientes com Leucemia

Dados

Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 23, página 182). Dados disponíveis em: <https://www.ime.usp.br/~giapaula/textoregressao.htm>, arquivo sobrev.dat.

O conjunto de dados a seguir refere-se à classificação de pacientes com leucemia segundo a ausência ou presença de uma característica morfológica nas células brancas. Os pacientes classificados de AG positivo foram aqueles identificados com a presença dessa característica e os pacientes classificados em AG negativo os que não apresentaram. O conjunto de dados também inclui o tempo de sobrevivência do paciente (em semanas) após o diagnóstico da doença e o número de células brancas (WBC) no momento do diagnóstico.

Variáveis:

• WBC: número de células brancas no momento do diagnóstico;

• TEMPO: tempo de sobrevivência do paciente (em semanas) após o diagnóstico da doença;

• AG: 0 em caso da ausência da característica e 1 em caso da presença.

Supondo que o tempo de sobrevivência após o diagnóstico segue uma distribuição gama, vamos propor um modelo para explicar o tempo médio de sobrevivência dados \(log(WBC)\) e AG (= 1 se positivo, = 0 se negativo).

Denotaremos por \(T_{ij}\) o tempo de sobrevivência do i-ésimo paciente com a j-ésima classificação, \(i = 1,…,33\) e \(j = 1,2\).

# Obtendo os dados
dados = read.csv("dados/sobrev.csv")
dados$AG = factor(dados$AG)
attach(dados)

Análise Descritiva

Médias estimadas para o tempo de sobrevivência, considerando o fator AG:

# AG Negativo
dados_ag0 = dados %>% filter(AG==0)
mean(dados_ag0$TEMPO)

## [1] 17.9375

# AG Positivo
dados_ag1 = dados %>% filter(AG==1)
mean(dados_ag1$TEMPO)

## [1] 62.47059

A partir da diferença entre as médias amostrais acima, podemos supor que os pacientes com o fator AG positivo sobreviveram por mais tempo que aqueles com fator AG negativo.

Coeficientes de variação do tempo de sobrevivência, considerando o fator AG:

# Tempo, AG = 0
sd(dados_ag0$TEMPO)/mean(dados_ag0$TEMPO) * 100

## [1] 113.1853

# Tempo, AG = 1
sd(dados_ag1$TEMPO)/mean(dados_ag1$TEMPO) * 100

## [1] 87.00598

Box-plots

Ignorando o fator AG, temos o seguinte boxplot:

plot_ly(type = "box", y = TEMPO, name = "TEMPO", color = I("slategray"))

Considerando o fator AG, temos os seguintes boxplots:

plot_ly(type = "box",
        y = TEMPO, x = AG, name = AG,
        color = AG, colors = c("darkred","darkgreen")) %>% 
  layout(title = "TEMPO x AG")

Os gráficos de box-plot evidenciam a assimetria na distribuição da variável tempo de sobrevivência e também apontam a presença de um outlier entre os indivíduos com fator AG negativo. Esse outlier corresponde ao paciente que, mesmo com o fator AG negativo, conseguiu sobreviver a 65 semanas. Nos dados, esse paciente é a 19ª observação.

Densidades, Dispersões e Correlações

Ignorando o fator AG, temos a seguinte densidade aproximada do tempo de sobrevivência:

ggplot(dados,aes(x=TEMPO)) +
  geom_density(fill = "slategray", color ="slategray",lwd = 1, alpha = 0.5) +
  ggtitle("Densidade: TEMPO") +
  theme_bw()

Considerando o fator AG, temos as seguintes dispersões, densidades e correlações:

g = ggpairs(dados[c("WBC","TEMPO")],
            aes(color = AG),
            lower = list(continuous = wrap("smooth",col="black")),
            diag = list(continuous = wrap("densityDiag", alpha = 0.5, size = 1)))
ggplotly(g)

De maneira geral, percebemos uma correlação fraca de -0,33 entre as variáveis WBC e TEMPO. Considerando o fator AG, a correlação é praticamente nula no caso de AG negativo e moderada no caso de AG positivo. Em outras palavras, na presença da característica morfológica, o tempo de sobrevivência tende a diminuir conforme o número de celulas brancas aumentam.

Ajuste

Supondo que \(T_{ij}\) são variáveis aleatórias independentes tais que \(T_{ij} \sim Gama(\mu_{ij},\phi)\) e que \(g(\mu_{ij}) = \eta_{ij}\), com \(\eta_{ij} = \bf{x}_{ij}^\top\beta\), \(\bf{x}_{ij}^\top\) contendo os valores das variáveis explicativas \(log(WBC)\) e \(AG\), e \(\beta\) o vetor de parâmetros.

Assumimos que \(T_{ij}\) segue uma distribuição Gama de média \(\mu_i\) e parâmetro de dispersão \(\phi^{-1}\). Vamos, então, propor um modelo cuja parte sistemática é dada por

\[ \log(\mu_{ij}) = \alpha + \gamma_j + \beta\log(WBC_{ij}) \; , \]

em que \(\gamma_{1} = 0\).

Modelo com Ligação Logarítmica

model1 = glm(TEMPO ~ log(WBC) + AG, dados, family = Gamma(link = "log"))
resumo = summary(model1)
resumo

## 
## Call:
## glm(formula = TEMPO ~ log(WBC) + AG, family = Gamma(link = "log"), 
##     data = dados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1746  -1.2663  -0.4251   0.4962   1.9048  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.8154     1.3487   4.312 0.000161 ***
## log(WBC)     -0.3044     0.1375  -2.213 0.034617 *  
## AG1           1.0176     0.3642   2.794 0.008982 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 1.087715)
## 
##     Null deviance: 58.138  on 32  degrees of freedom
## Residual deviance: 40.319  on 30  degrees of freedom
## AIC: 301.49
## 
## Number of Fisher Scoring iterations: 8

Para o modelo obtido, ao nível de 5% de significância, as variáveis log(WBC) e AG são significativas.

shape = MASS::gamma.shape(model1)

desvio_escal = resumo$deviance * shape$alpha
nivel_descrit = 1 - pchisq(desvio_escal, model1$df.residual)
nivel_descrit

## [1] 0.1415408

Também ao nível de 5% de significância, não rejeitamos a hipótese de que o modelo seja adequado.

Seleção de Variáveis

MASS::stepAIC(model1)

## Start:  AIC=301.49
## TEMPO ~ log(WBC) + AG
## 
##            Df Deviance    AIC
## <none>          40.319 301.49
## - log(WBC)  1   46.198 304.89
## - AG        1   47.808 306.37

## 
## Call:  glm(formula = TEMPO ~ log(WBC) + AG, family = Gamma(link = "log"), 
##     data = dados)
## 
## Coefficients:
## (Intercept)     log(WBC)          AG1  
##      5.8154      -0.3044       1.0176  
## 
## Degrees of Freedom: 32 Total (i.e. Null);  30 Residual
## Null Deviance:       58.14 
## Residual Deviance: 40.32     AIC: 301.5

Pelo critério de seleção de Akaike, devemos manter as duas variáveis no modelo.

Diagnóstico

Envelope

# ENVELOPE
fit.model = model1
source("envel_gama")

Pelo gráfico de envelope, verificamos que todos os pontos se encontram dentro das bandas.

par(mfrow=c(2,2))

# RESÍDUOS STUDENTIZADOS
residualPlot(model1, type = "deviance", id = TRUE,
             ylab = "res. compon. do desvio", xlab = "valor ajustado",
             pch = 16, smooth = F)
abline(h = c(-2,2), col="red")

# VARIÁVEL Z
w = model1$weights
eta = model1$linear.predictors
z = eta + resid(model1, type="pearson")/sqrt(w)
plot(model1$linear.predictors, z,
     xlab = "Preditor Linear", ylab = "Variável z", pch = 16)

# ALAVANCAGEM
cut = min(sort(hatvalues(model1), decreasing = TRUE)[1:3])
plot(hatvalues(model1), pch = 16,
     xlab = "observação", ylab = "medida h", main = "Alavancagem")
text(hatvalues(model1) * as.integer(hatvalues(model1) >= cut),
     cex=0.8, p=1, offset=0.3)

# DISTÂNCIA DE COOK
plot(model1, which=4, lwd=3, main = "Distância de Cook", caption = F)

No gráfico de resíduos por valores ajustados, percebemos uma boa distribuição dos pontos no intervalo de -2 a 2. Embora seja perceptível também um afunilamento à direita, a quantidade de pontos nessa região não é suficiente para indicar que seja uma tendência. Não há presença de pontos aberrantes, uma vez que todos os pontos estão dentro ou muito próximo do intervalo de -2 a 2. A observação 33 fica em evidência no gráfico de distância de Cook, mas não a consideraremos para análise uma vez que a distância não é expressiva em escala. No gráfico de alavancagem, a observação 2 se destaca e, sendo assim, investigaremos sua possível influência na qualidade do ajuste.

Pontos Influentes

summary(influence.measures(model1))

## Potentially influential observations of
##   glm(formula = TEMPO ~ log(WBC) + AG, family = Gamma(link = "log"),      data = dados) :
## 
##   dfb.1_ dfb.l(WB dfb.AG1 dffit cov.r   cook.d hat  
## 2  0.09  -0.09     0.04    0.11  1.37_*  0.01   0.20

As medidas de influência acima apontam apenas para a observação 2 como possível ponto influente.

Impacto das remoções das observações

model1_sem2 = glm(TEMPO ~ log(WBC) + AG, dados, subset = -2,
                   family = Gamma(link = "log"))
impacto(model1, model1_sem2)

##             coef_com_i coef_sem_i   impacto   pval_com_i  pval_sem_i
## (Intercept)  5.8154109  5.6913329 2.1336073 0.0001610584 0.000584179
## log(WBC)    -0.3044004 -0.2919183 4.1005537 0.0346169527 0.062619407
## AG1          1.0176452  1.0090587 0.8437636 0.0089824909 0.011413245

Ao nível de 5% de significância, verificamos que a variável log(WBC) deixa de ser significativa quando removemos a observação 2. Sendo assim, ela causa mudança infefrencial e classifica-se como ponto de alavanca e influente.

Interpretação

Com essa análise, chegamos ao seguinte modelo para explicar o tempo médio de sobrevivência:

\[ \begin{aligned} \log(\mu_{i1}) &= 5.8 - 0.3\log(WBC_{i1}) \; ; \\ \log(\mu_{i2}) &= 6.8 - 0.3\log(WBC_{i2}) \; . \end{aligned} \]

Fixando-se o fator \(\gamma_{j}\), temos que

\[ \frac{\mu_{j}(x+1)}{\mu_{j}(x)} = \frac{\exp{(5.8+ \gamma_j -0.3(x+1))}}{\exp{(5.8 + \gamma_j - 0.3x)}} =\exp{(-0.3) = 0.74} \; . \]

Alterando-se o fator \(\gamma_j\), temos que

\[ \frac{\mu_{i2}}{\mu_{i1}} = \frac{\exp{(5.8 + \gamma_2 - 0.3X_{i2})}}{\exp{(5.8 + \gamma_1 - 0.3X_{i1})}} = \exp{(\gamma_2)} = 2.72 \]

De onde concluímos que, ao fixarmos o fator característica morfológica, o acréscimo de 1 unidade à variável \(\log(WBC)\) representa, em média, uma redução de 26% do tempo de sobrevivência. Já considerando a característica morfológica, o fato do paciente apresentá-la representa um aumento médio de 172% no tempo de sobrevivência.

III. Modelo para Dados Binários

1. Dose-Resposta

Dados

Os conjuntos de dados apresentados a seguir são provenientes de um experimento dose-resposta, conduzido para avaliar a influência dos extratos vegetais “aquoso frio de folhas”, “aquoso frio de frutos” e um extrato químico, respectivamente, na morte de um determinado tipo de caramujo. Para cada conjunto, ajustaremos um modelo logístico linear simples e um modelo log-log linear simples. Para o melhor ajuste, usando envelope como critério, encontraremos um intervalo assintótico de 95% para a dose letal \(DL_{50}\).

Dados do Extrato Aquoso Frio de Folhas:

dose1 = read.table("dados/dose1.dat", col.names = c("dose","expostos","mortos"))
dose1 = dose1 %>% mutate(restantes = expostos - mortos,
                         prop.mort = mortos/expostos,
                         prop.rest = restantes/expostos)
dose1

##   dose expostos mortos restantes prop.mort prop.rest
## 1    0       50      4        46      0.08      0.92
## 2   15       50      5        45      0.10      0.90
## 3   20       50     14        36      0.28      0.72
## 4   25       50     29        21      0.58      0.42
## 5   30       50     38        12      0.76      0.24
## 6   35       50     41         9      0.82      0.18
## 7   40       50     47         3      0.94      0.06

Dados do Extrato Aquoso Frio de Frutos:

dose2 = read.table("dados/dose2.dat", col.names = c("dose","expostos","mortos"))
dose2 = dose2 %>% mutate(restantes = expostos - mortos,
                         prop.mort = mortos/expostos,
                         prop.rest = restantes/expostos)
dose2

##   dose expostos mortos restantes  prop.mort  prop.rest
## 1    0       65      2        63 0.03076923 0.96923077
## 2  100       52      2        50 0.03846154 0.96153846
## 3  150       52      8        44 0.15384615 0.84615385
## 4  200       52     20        32 0.38461538 0.61538462
## 5  250       52     41        11 0.78846154 0.21153846
## 6  300       52     48         4 0.92307692 0.07692308
## 7  350       52     51         1 0.98076923 0.01923077
## 8  400       52     52         0 1.00000000 0.00000000

Dados do Extrato Químico:

dose3 = read.table("dados/dose3.dat", col.names = c("dose","expostos","mortos"))
dose3 = dose3 %>% mutate(restantes = expostos - mortos,
                         prop.mort = mortos/expostos,
                         prop.rest = restantes/expostos)
dose3

##    dose expostos mortos restantes  prop.mort prop.rest
## 1 0.000       30      2        28 0.06666667 0.9333333
## 2 0.010       30      4        26 0.13333333 0.8666667
## 3 0.015       30      4        26 0.13333333 0.8666667
## 4 0.020       30     11        19 0.36666667 0.6333333
## 5 0.030       29      9        20 0.31034483 0.6896552
## 6 0.040       30     14        16 0.46666667 0.5333333
## 7 0.050       30     24         6 0.80000000 0.2000000

1.1 Proporção de Caramujos Mortos por Dose

p1 = ggplot(dose1, aes(x = dose, y = prop.mort)) +
  geom_col(fill = "darkred") +
  ylab("prop. mortos") +
  labs(title = "Extrato Aquoso Frio de Folhas") +
  theme_bw()

p2 = ggplot(dose2, aes(x = dose, y = prop.mort)) +
  geom_col(fill = "darkred") +
  ylab("prop. mortos") +
  labs(title = "Extrato Aquoso Frio de Frutos") +
  theme_bw()

p3 = ggplot(dose3, aes(x = dose, y = prop.mort)) +
  geom_col(fill = "darkred") +
  ylab("prop. mortos") +
  labs(title = "Extrato Químico") +
  theme_bw()

plot_grid(p1,p2,p3,ncol=2)

1.2 Modelos Logístico Linear Simples

# EAF Folhas
log1 = glm(as.matrix(dose1[3:4]) ~ dose, data = dose1, family=binomial)
summary(log1)

## 
## Call:
## glm(formula = as.matrix(dose1[3:4]) ~ dose, family = binomial, 
##     data = dose1)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7  
##  2.1850  -1.7488  -0.9036   0.7107   0.7624  -0.4650   0.4818  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.80668    0.45190  -8.424   <2e-16 ***
## dose         0.15707    0.01707   9.199   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 172.446  on 6  degrees of freedom
## Residual deviance:  10.184  on 5  degrees of freedom
## AIC: 40.101
## 
## Number of Fisher Scoring iterations: 5

# EAF Frutos
log2 = glm(as.matrix(dose2[3:4]) ~ dose, data = dose2, family=binomial)
summary(log2)

## 
## Call:
## glm(formula = as.matrix(dose2[3:4]) ~ dose, family = binomial, 
##     data = dose2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8190  -0.3956   0.1640   0.7142   2.1973  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.546200   0.565304  -9.811   <2e-16 ***
## dose         0.026539   0.002561  10.363   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 353.3692  on 7  degrees of freedom
## Residual deviance:   7.0361  on 6  degrees of freedom
## AIC: 33.504
## 
## Number of Fisher Scoring iterations: 5

# E Químico
log3 = glm(as.matrix(dose3[3:4]) ~ dose, data = dose3, family=binomial)
summary(log3)

## 
## Call:
## glm(formula = as.matrix(dose3[3:4]) ~ dose, family = binomial, 
##     data = dose3)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7  
## -0.05150   0.03496  -0.65291   1.61394  -0.83739  -1.02803   0.98908  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.6016     0.3662  -7.105  1.2e-12 ***
## dose         71.0948    11.3952   6.239  4.4e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.0428  on 6  degrees of freedom
## Residual deviance:  5.7713  on 5  degrees of freedom
## AIC: 33.346
## 
## Number of Fisher Scoring iterations: 4

Testes de Ajuste

# EAF Folhas
1 - pchisq(log1$deviance, log1$df.residual)

## [1] 0.07019215

# EAF Frutos
1 - pchisq(log2$deviance, log1$df.residual)

## [1] 0.2179659

# E Químico
1 - pchisq(log3$deviance, log1$df.residual)

## [1] 0.3291096

1.3 Modelos Complementar Log-log

# EAF Folhas
cloglog1 = glm(as.matrix(dose1[3:4]) ~ dose, data = dose1,
               family = binomial(link = "cloglog"))
summary(cloglog1)

## 
## Call:
## glm(formula = as.matrix(dose1[3:4]) ~ dose, family = binomial(link = "cloglog"), 
##     data = dose1)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7  
##  0.9330  -2.0977  -0.7403   1.3129   1.3331  -0.5276  -0.5728  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.99021    0.30754  -9.723   <2e-16 ***
## dose         0.10348    0.01017  10.175   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 172.4463  on 6  degrees of freedom
## Residual deviance:   9.9264  on 5  degrees of freedom
## AIC: 39.844
## 
## Number of Fisher Scoring iterations: 4

# EAF Frutos
cloglog2 = glm(as.matrix(dose2[3:4]) ~ dose, data = dose2,
               family = binomial(link = "cloglog"))
summary(cloglog2)

## 
## Call:
## glm(formula = as.matrix(dose2[3:4]) ~ dose, family = binomial(link = "cloglog"), 
##     data = dose2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.48238  -0.94709  -0.03929   0.18022   1.80452  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.988340   0.372324  -10.71   <2e-16 ***
## dose         0.016428   0.001498   10.96   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 353.3692  on 7  degrees of freedom
## Residual deviance:   8.6014  on 6  degrees of freedom
## AIC: 35.069
## 
## Number of Fisher Scoring iterations: 7

# E Químico
cloglog3 = glm(as.matrix(dose3[3:4]) ~ dose, data = dose3,
               family = binomial(link = "cloglog"))
summary(cloglog3)

## 
## Call:
## glm(formula = as.matrix(dose3[3:4]) ~ dose, family = binomial(link = "cloglog"), 
##     data = dose3)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7  
## -0.25501  -0.02679  -0.62490   1.74222  -0.60022  -0.90515   0.58896  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.4987     0.3128  -7.988 1.38e-15 ***
## dose         56.7822     8.4721   6.702 2.05e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.043  on 6  degrees of freedom
## Residual deviance:  5.018  on 5  degrees of freedom
## AIC: 32.593
## 
## Number of Fisher Scoring iterations: 4

Teste de Ajuste

# EAF Folhas
1 - pchisq(cloglog1$deviance, cloglog1$df.residual)

## [1] 0.07734658

# EAF Frutos
1 - pchisq(cloglog2$deviance, cloglog2$df.residual)

## [1] 0.1972691

# E químico
1 - pchisq(cloglog3$deviance, cloglog3$df.residual)

## [1] 0.4136878

1.4 Envelopes

Extrato Aquoso Frio de Folhas:

fit.model = log1
ntot = dose1$expostos
source("envelr_bino.txt")

fit.model = cloglog1
ntot = dose1$expostos
source("envelr_bino_cloglog.txt")

Extrato Aquoso Frio de Frutos:

fit.model = log2
ntot = dose2$expostos
source("envelr_bino.txt")

fit.model = cloglog2
ntot = dose2$expostos
source("envelr_bino_cloglog.txt")

Extrato Químico:

fit.model = log3
ntot = dose3$expostos
source("envelr_bino.txt")

fit.model = cloglog3
ntot = dose3$expostos
source("envelr_bino_cloglog.txt")

1.5 Curvas Logísticas Ajustadas

pi_log1 = exp(log1$coeff[1] + log1$coeff[2] * dose1$dose) /
  (1 + exp(log1$coeff[1] + log1$coeff[2] * dose1$dose))

p_log1 = ggplot(dose1, aes(x = dose, y = prop.mort)) +
  geom_point() +
  geom_line(y = pi_log1) +
  ylab("prop. mortos") +
  ggtitle("E.A.F.Folhas - Log simples") +
  theme_bw()

pi_cloglog1 = exp(cloglog1$coeff[1] + cloglog1$coeff[2]*dose1$dose) /
  (1 + exp(cloglog1$coeff[1] + cloglog1$coeff[2]*dose1$dose))

p_cloglog1 = ggplot(dose1, aes(x = dose, y = prop.mort)) +
  geom_point() +
  geom_line(y = pi_cloglog1) +
  ylab("prop. mortos") +
  ggtitle("E.A.F.Folhas - Complem. Log-log") +
  theme_bw()



pi_log2 = exp(log2$coeff[1] + log2$coeff[2] * dose2$dose) /
  (1 + exp(log2$coeff[1] + log2$coeff[2] * dose2$dose))

p_log2 = ggplot(dose2, aes(x = dose, y = prop.mort)) +
  geom_point() +
  geom_line(y = pi_log2) +
  ylab("prop. mortos") +
  ggtitle("E.A.F.Frutos - Log simples") +
  theme_bw()

pi_cloglog2 = exp(cloglog2$coeff[1] + cloglog2$coeff[2]*dose2$dose) /
  (1 + exp(cloglog2$coeff[1] + cloglog2$coeff[2]*dose2$dose))

p_cloglog2 = ggplot(dose2, aes(x = dose, y = prop.mort)) +
  geom_point() +
  geom_line(y = pi_cloglog2) +
  ylab("prop. mortos") +
  ggtitle("E.A.F.Frutos - Complem. Log-log") +
  theme_bw()



pi_log3 = exp(log3$coeff[1] + log3$coeff[2] * dose3$dose) /
  (1 + exp(log3$coeff[1] + log3$coeff[2] * dose3$dose))

p_log3 = ggplot(dose3, aes(x = dose, y = prop.mort)) +
  geom_point() +
  geom_line(y = pi_log3) +
  ylab("prop. mortos") +
  ggtitle("E. Químico - Log simples") +
  theme_bw()

pi_cloglog3 = exp(cloglog3$coeff[1] + cloglog3$coeff[2]*dose3$dose) /
  (1 + exp(cloglog3$coeff[1] + cloglog3$coeff[2]*dose3$dose))

p_cloglog3 = ggplot(dose3, aes(x = dose, y = prop.mort)) +
  geom_point() +
  geom_line(y = pi_cloglog3) +
  ylab("prop. mortos") +
  ggtitle("E. Químico - Complem. Log-log") +
  theme_bw()


plot_grid(p_log1,p_cloglog1,
          p_log2,p_cloglog2,
          p_log3,p_cloglog3, ncol = 2)

1.6 Dose Letal (\(DL_{50}\)) para os Modelos Log Linear Simples

# DOSE LETAL ESTIMATIVA PONTUAL
DL_100p = function(model,p){
  (1/model$coeff[2])*(log(p/(1-p)) - model$coeff[1])
}


# DOSE LETAL ESTIMATIVA INTERVALAR, 95% DE CONFIANÇA
DL_IC_log = function(model, p, gama){
  db = cbind(-1/model$coeff[2],
             1/(model$coeff[2]^2)*(model$coeff[1]-log(p/(1-p))))
  
  X = model.matrix(model)
  w = model$weights
  W = diag(w)
  
  Var = db %*% solve(t(X) %*% W %*% X) %*% t(db)
  
  LI = DL_100p(model, 0.5) - qnorm((1+gama)/2)*sqrt(Var)
  LS = DL_100p(model, 0.5) + qnorm((1+gama)/2)*sqrt(Var)
  
  return(c(LI,LS))
}

data.frame(extrato = c("EAF de Folhas","EAF de Frutos", "E Químico"),
           LI = c(DL_IC_log(log1, 0.5, 0.95)[1],
                  DL_IC_log(log2, 0.5, 0.95)[1],
                  DL_IC_log(log3, 0.5, 0.95)[1]),
           DL_100p = c(DL_100p(log1,0.5),DL_100p(log2,0.5),DL_100p(log3,0.5)),
           LS = c(DL_IC_log(log1, 0.5, 0.95)[2],
                  DL_IC_log(log2, 0.5, 0.95)[2],
                  DL_IC_log(log3, 0.5, 0.95)[2]))

##         extrato           LI      DL_100p           LS
## 1 EAF de Folhas  22.50628769  24.23566646  25.96504522
## 2 EAF de Frutos 197.03961411 208.98193852 220.92426293
## 3     E Químico   0.03131498   0.03659294   0.04187089

2. Eficiência do Tratamento de Pacientes com Leucemia

2.1 Dados

Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 20, página 277). Dados disponíveis em: <https://www.ime.usp.br/~giapaula/textoregressao.htm>, arquivo leuc.dat.

Os dados a seguir referem-se a um estudo com 51 pacientes adultos, previamente diagnosticados com um tipo agudo de leucemia, e que receberam um tipo de tratamento. A eficiência ou não do tratamento foi verificada após um certo período. As variáveis em estudo são:

• idade (\(x_1\)): idade do paciente (em anos);

• mancha (\(x_2\)): mancha diferencial da doença (em %);

• infilt (\(x_3\)): infiltração na medula (em %);

• leuc (\(x_4\)): células com leucemia (em %);

• malig (\(x_5\)): malignidade da doença (\(\times 10^3\));

• tmax (\(x_6\)): temperatura máxima pré-tratamento (\(\times 10^\circ F\));

• trat (\(y\)): tratamento (0: não-satisfatório, 1: satisfatório).

Consideremos o seguinte modelo logístico linear, para explicar a probabilidade de eficiência do tratamento a partir das seis primeiras variáveis explicativas acima:

\[ Y \sim Bernoulli(\pi(\mathbf{x})) \]

\[ \log \left\{ \frac{\pi(\mathbf{x})}{1-\pi(\mathbf{x})} \right\} = \beta_1 + \beta_2x_2 + ... + \beta_6x_6 \]

dados = read.table("dados/leuce.dat",
                   col.names = c("idade","mancha","infilt",
                                 "leuc","malig","tmax",
                                 "trat","tvida","status"))
dados$trat = factor(dados$trat)
dados = dados[1:7]
attach(dados)

Box-plots

p1 = dados %>% ggplot(aes(x = trat, y = idade, fill = trat)) +
  geom_boxplot() +
  scale_fill_manual(values=c("#BC3C29", "#20854E")) +
  ggtitle("idade x tratamento") +
  theme(legend.position = "none")

p2 = dados %>% ggplot(aes(x = trat, y = mancha, fill = trat)) +
  geom_boxplot() +
  scale_fill_manual(values=c("#BC3C29", "#20854E")) +
  ggtitle("mancha diferencial x tratamento") +
  theme(legend.position = "none")

p3 = dados %>% ggplot(aes(x = trat, y = infilt, fill = trat)) +
  geom_boxplot() +
  scale_fill_manual(values=c("#BC3C29", "#20854E")) +
  ggtitle("infiltração x tratamento") +
  theme(legend.position = "none")

p4 = dados %>% ggplot(aes(x = trat, y = leuc, fill = trat)) +
  geom_boxplot() +
  scale_fill_manual(values=c("#BC3C29", "#20854E")) +
  ggtitle("células leucêmicas x tratamento") +
  theme(legend.position = "none")

p5 = dados %>% ggplot(aes(x = trat, y = malig, fill = trat)) +
  geom_boxplot() +
  scale_fill_manual(values=c("#BC3C29", "#20854E")) +
  ggtitle("malignidade x tratamento") +
  theme(legend.position = "none")

p6 = dados %>% ggplot(aes(x = trat, y = tmax, fill = trat)) +
  geom_boxplot() +
  scale_fill_manual(values=c("#BC3C29", "#20854E")) +
  ggtitle("temperatura máxima x tratamento") +
  theme(legend.position = "none")

grid.arrange(p1,p2,p3,p4,p5,p6)

Observando os box-plots acima, fica evidente que o tratamento foi mais eficaz em pacientes mais jovens e naqueles com maior porcentagem de células leucêmicas.

Correlograma

corrplot::corrplot(cor(dados[1:6]),  method = "color", type = "lower",
                   addCoef.col = "white", tl.col = "black")

Dentre as variáveis explicativas, o par mancha diferencial e infiltração na medula foi o único com alto grau de correlação. Então, vale estudar modelos considerando as duas individualmente junto as demais.

2.2 Modelo Logístico Linear

Modelo Completo

model = glm(trat ~ idade + mancha + infilt + leuc + malig + tmax,
            family = binomial)
summary(model)

## 
## Call:
## glm(formula = trat ~ idade + mancha + infilt + leuc + malig + 
##     tmax, family = binomial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.56756  -0.55762  -0.05269   0.59038   2.38751  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 108.33115   41.84379   2.589  0.00963 **
## idade        -0.06231    0.02746  -2.269  0.02327 * 
## mancha       -0.00469    0.04005  -0.117  0.90677   
## infilt        0.03104    0.03789   0.819  0.41264   
## leuc          0.37281    0.13247   2.814  0.00489 **
## malig         0.03267    0.04605   0.710  0.47801   
## tmax         -0.11162    0.04263  -2.618  0.00884 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 70.524  on 50  degrees of freedom
## Residual deviance: 39.275  on 44  degrees of freedom
## AIC: 53.275
## 
## Number of Fisher Scoring iterations: 6

Para o modelo com todas as variáveis fornecidas, mancha, infiltração e malignidade não foram significativas.

Seleção de Variáveis

MASS::stepAIC(model)

## Start:  AIC=53.28
## trat ~ idade + mancha + infilt + leuc + malig + tmax
## 
##          Df Deviance    AIC
## - mancha  1   39.289 51.289
## - infilt  1   40.014 52.014
## - malig   1   40.115 52.115
## <none>        39.275 53.275
## - idade   1   45.737 57.737
## - tmax    1   48.760 60.760
## - leuc    1   52.213 64.213
## 
## Step:  AIC=51.29
## trat ~ idade + infilt + leuc + malig + tmax
## 
##          Df Deviance    AIC
## - malig   1   40.136 50.136
## - infilt  1   41.042 51.042
## <none>        39.289 51.289
## - idade   1   45.784 55.784
## - tmax    1   48.882 58.882
## - leuc    1   52.728 62.728
## 
## Step:  AIC=50.14
## trat ~ idade + infilt + leuc + tmax
## 
##          Df Deviance    AIC
## <none>        40.136 50.136
## - infilt  1   43.265 51.265
## - idade   1   46.438 54.438
## - tmax    1   48.971 56.971
## - leuc    1   57.602 65.602

## 
## Call:  glm(formula = trat ~ idade + infilt + leuc + tmax, family = binomial)
## 
## Coefficients:
## (Intercept)        idade       infilt         leuc         tmax  
##    95.56766     -0.06026      0.03413      0.40673     -0.09944  
## 
## Degrees of Freedom: 50 Total (i.e. Null);  46 Residual
## Null Deviance:       70.52 
## Residual Deviance: 40.14     AIC: 50.14

No método stepwise AIC, as variáveis selecionadas foram idade, infilt, leuc e tmax. Como infilt e mancha são fortemente correlacionadas, vamos testar também um modelo com mancha ao invés de infilt.

Modelos Reduzidos

# Com a variável INFILT
modelr1 = glm(trat ~ idade + infilt + leuc + tmax, family = binomial)
summary(modelr1)

## 
## Call:
## glm(formula = trat ~ idade + infilt + leuc + tmax, family = binomial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.73886  -0.56473  -0.05442   0.62185   2.26516  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 95.56766   38.59482   2.476  0.01328 * 
## idade       -0.06026    0.02678  -2.250  0.02445 * 
## infilt       0.03413    0.02079   1.641  0.10077   
## leuc         0.40673    0.13034   3.121  0.00181 **
## tmax        -0.09944    0.03954  -2.515  0.01191 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 70.524  on 50  degrees of freedom
## Residual deviance: 40.136  on 46  degrees of freedom
## AIC: 50.136
## 
## Number of Fisher Scoring iterations: 6

# Com a variável MANCHA
modelr2 = glm(trat ~ idade + mancha + leuc + tmax, family = binomial)
summary(modelr2)

## 
## Call:
## glm(formula = trat ~ idade + mancha + leuc + tmax, family = binomial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.87043  -0.58810  -0.05384   0.65397   2.30025  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 94.45927   37.25941   2.535  0.01124 * 
## idade       -0.05637    0.02632  -2.141  0.03225 * 
## mancha       0.03029    0.02191   1.382  0.16682   
## leuc         0.41243    0.12938   3.188  0.00143 **
## tmax        -0.09854    0.03826  -2.576  0.01000 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 70.524  on 50  degrees of freedom
## Residual deviance: 41.222  on 46  degrees of freedom
## AIC: 51.222
## 
## Number of Fisher Scoring iterations: 6

Os modelo reduzido com a variável infilt quanto o com a variável mancha apresentam as mesmas estimativas aproximadas para os parâmetros. Também nos dois, essas variáveis permanecem não significativas. Adiante, estudaremos uma possível influência de pontos sobre a inferência do modelo e avaliaremos a remoção das variáveis.

Testes de Ajuste

1 - pchisq(modelr1$deviance, modelr1$df.residual)

## [1] 0.7153545

1 - pchisq(modelr2$deviance, modelr2$df.residual)

## [1] 0.6723047

Ao nível de 5% de significância, não rejeitamos a hipótese de que os modelos sejam adequados.

Envelopes

fit.model = modelr1
source("envel_bino.txt")

fit.model = modelr2
source("envel_bino.txt")

Pelos gráficos de envelope, temos que o modelo reduzido com a variável mancha ficou melhor ajustado e, então, o escolheremos para prosseguir nossa análise.

2.3 Diagnóstico

par(mfrow=c(2,2))

# ALAVANCAGEM
cut = min(sort(hatvalues(modelr2), decreasing = TRUE)[1:3])
plot(hatvalues(modelr2), pch = 16,
     xlab = "observação", ylab = "medida h", main = "Alavancagem")
text(hatvalues(modelr2) * as.integer(hatvalues(modelr2) >= cut),
     cex=0.8, p=1, offset=0.4)

# DISTÂNCIA DE COOK
plot(modelr2, which=4, lwd=3, main = "Distância de Cook", caption = F)

# RESÍDUO COMPONENTE DO DESVIO POR ÍNDICE
plot(residuals(modelr2), pch = 16,
     ylab = "res. compon. do desvio", xlab = "índice")
abline(h = c(-2,2), col="red")

No gráfico de resíduo componente do desvio, verificamos que os pontos estão bem distribuídos aleatoriamente e sem nenhuma observação aberrante muito além dos limites de confiança. Nos gráficos de alavancagem e distância de Cook, destacam–se os pontos 2, 13, 32, 47 e 48. A seguir, verificaremos se há influência desses pontos sobre a inferência do modelo.

Pontos Influentes

summary(influence.measures(modelr2))

## Potentially influential observations of
##   glm(formula = trat ~ idade + mancha + leuc + tmax, family = binomial) :
## 
##    dfb.1_ dfb.idad dfb.mnch dfb.leuc dfb.tmax dffit   cov.r   cook.d hat    
## 2  -0.61  -0.36    -0.22     0.02     0.62     1.06_*  1.31    0.14   0.33_*
## 13 -0.09   0.23     0.60     0.14     0.05    -0.74    1.33_*  0.06   0.28  
## 32  0.23  -0.16    -0.44    -0.46    -0.19    -1.03_*  0.87    0.20   0.20  
## 47 -0.45   0.58    -0.22    -0.40     0.44     0.71    0.57_*  0.20   0.07  
## 48 -0.30  -0.34     0.55    -0.44     0.29    -1.11_*  1.11    0.18   0.28

Impacto dos Pontos Influentes

Observação 2

modelr2_sem2 = glm(trat ~ idade + mancha + leuc + tmax, subset = -2,
                   family = binomial)
impacto(modelr2, modelr2_sem2)[3:5]

##                 impacto  pval_com_i  pval_sem_i
## (Intercept) -19.7958316 0.011238999 0.007759372
## idade        12.2404843 0.032251493 0.061370014
## mancha      -13.4866607 0.166824037 0.127772519
## leuc         -0.7081507 0.001434345 0.001634411
## tmax        -19.8367890 0.010004635 0.007094660

A 5% de significância, removendo-se a observação 2, a variável idade deixa de ser significativa.

Observação 13

modelr2_sem13 = glm(trat ~ idade + mancha + leuc + tmax, subset = -13,
                   family = binomial)
impacto(modelr2, modelr2_sem13)[3:5]

##               impacto  pval_com_i  pval_sem_i
## (Intercept) -3.204997 0.011238999 0.009617055
## idade       -8.835451 0.032251493 0.024988565
## mancha      32.179710 0.166824037 0.393890235
## leuc         2.477318 0.001434345 0.001728282
## tmax        -1.989041 0.010004635 0.008939169

A 5% de significância, removendo-se a observação 13, não há mudança inferencial no modelo.

Observação 32

modelr2_sem32 = glm(trat ~ idade + mancha + leuc + tmax, subset = -32,
                   family = binomial)
impacto(modelr2, modelr2_sem32)[3:5]

##                impacto  pval_com_i  pval_sem_i
## (Intercept)   8.654518 0.011238999 0.026129747
## idade         7.325326 0.032251493 0.052668960
## mancha      -36.997935 0.166824037 0.091633885
## leuc        -17.602426 0.001434345 0.001022189
## tmax          6.768428 0.010004635 0.021156787

A 5% de significância, removendo-se a observação 32, a variável idade deixa de ser significativa.

Observação 47

modelr2_sem47 = glm(trat ~ idade + mancha + leuc + tmax, subset = -47,
                   family = binomial)
impacto(modelr2, modelr2_sem47)[3:5]

##               impacto  pval_com_i  pval_sem_i
## (Intercept) -33.67437 0.011238999 0.006535694
## idade       -50.54577 0.032251493 0.013157635
## mancha      -29.65272 0.166824037 0.103656074
## leuc        -26.05928 0.001434345 0.001848562
## tmax        -33.01815 0.010004635 0.006022423

A 5% de significância, removendo-se a observação 47, não há mudança inferencial no modelo.

Observação 48

modelr2_sem48 = glm(trat ~ idade + mancha + leuc + tmax, subset = -48,
                   family = binomial)
impacto(modelr2, modelr2_sem48)[3:5]

##               impacto  pval_com_i  pval_sem_i
## (Intercept) -13.22389 0.011238999 0.009203440
## idade        12.27356 0.032251493 0.061894967
## mancha       33.26921 0.166824037 0.371617432
## leuc        -14.88983 0.001434345 0.001479636
## tmax        -12.88917 0.010004635 0.008270069

A 5% de significância, removendo-se a observação 48, a variável idade deixa de ser significativa.

2.4 Modelo Logístico Linear - Reajuste

As variáveis mancha e infilt não foram significativas em nenhuma circustância vista anteriormente. Então, podemos considerar removê-las. Além disso, a variável idade deixou de ser significativa com a remoção de algumas observações. Então, testaremos a seguir dois modelos sem as variáveis mancha e infilt: um com a variável idade e outro sem.

# Com a variável IDADE
modelr3 = glm(trat ~ idade + leuc + tmax, family = binomial)
summary(modelr3)

## 
## Call:
## glm(formula = trat ~ idade + leuc + tmax, family = binomial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.76104  -0.68683  -0.09747   0.67388   2.16510  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 87.38804   35.45816   2.465  0.01372 * 
## idade       -0.05850    0.02558  -2.287  0.02218 * 
## leuc         0.38493    0.12152   3.168  0.00154 **
## tmax        -0.08897    0.03607  -2.467  0.01363 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 70.524  on 50  degrees of freedom
## Residual deviance: 43.265  on 47  degrees of freedom
## AIC: 51.265
## 
## Number of Fisher Scoring iterations: 6

# Sem a variável IDADE
modelr4 = glm(trat ~ leuc + tmax, family = binomial)
summary(modelr4)

## 
## Call:
## glm(formula = trat ~ leuc + tmax, family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2508  -0.9029  -0.1332   0.7621   1.7328  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) 70.82281   30.70009   2.307  0.02106 * 
## leuc         0.33631    0.10422   3.227  0.00125 **
## tmax        -0.07465    0.03144  -2.374  0.01758 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 70.524  on 50  degrees of freedom
## Residual deviance: 49.645  on 48  degrees of freedom
## AIC: 55.645
## 
## Number of Fisher Scoring iterations: 5

Em ambos os modelos, a 5% de significâncias, todas as variáveis são significativas.

Testes de Ajuste

1 - pchisq(modelr3$deviance, modelr3$df.residual)

## [1] 0.6280164

1 - pchisq(modelr4$deviance, modelr4$df.residual)

## [1] 0.4075115

Ao nível de 5% de significância, não rejeitamos a hipótese de que os modelos sejam adequados.

Envelopes

fit.model = modelr3
source("envel_bino.txt")

fit.model = modelr4
source("envel_bino.txt")

Pelos gráficos de envelope, temos que os dois modelos ficaram bem ajustados. Então, pelo princípio da parcimônia, optaremos pelo modelo sem a variável idade.

Qualidade do Ajuste - Hosmer-Lemeshow

HosmerLemeshowTest(fit=fitted(modelr4), trat, ngr=10, X=cbind(leuc, tmax))$gof

## 
##  le Cessie-van Houwelingen-Copas-Hosmer global goodness of fit test
## 
## data:  fitted(modelr4) and trat
## z = 0.006054, p-value = 0.9952

Pelo teste de Homser-Lemeshow, ao nível de 5% de significância, não rejeitamos a hipótese de bom ajuste do modelo.

2.5 Diagnóstico

par(mfrow=c(2,2))

# ALAVANCAGEM
cut = min(sort(hatvalues(modelr4), decreasing = TRUE)[1:3])
plot(hatvalues(modelr4), pch = 16,
     xlab = "observação", ylab = "medida h", main = "Alavancagem")
text(hatvalues(modelr4) * as.integer(hatvalues(modelr4) >= cut),
     cex=0.8, p=1, offset=0.4)

# DISTÂNCIA DE COOK
plot(modelr4, which=4, lwd=3, main = "Distância de Cook", caption = F)

# RESÍDUO COMPONENTE DO DESVIO POR ÍNDICE
plot(residuals(modelr4), pch = 16,
     ylab = "res. compon. do desvio", xlab = "índice")
abline(h = c(-2,2), col="red")

Pelo gráfico de resíduo componente do desvio, verificamos que os pontos estão bem distribuídos aleatoriamente no intervalo (-2,2) e sem pontos aberrantes para muito além dele. Nos gráficos de alavancagem e distância de Cook, destacam-se os pontos 2, 32, 40 e 42. Então, a seguir vamos verificar a se há influência desses pontos sobre a inferência do modelo.

Pontos Influentes

summary(influence.measures(modelr4))

## Potentially influential observations of
##   glm(formula = trat ~ leuc + tmax, family = binomial) :
## 
##    dfb.1_ dfb.leuc dfb.tmax dffit   cov.r   cook.d hat    
## 2  -0.72   0.02     0.71     0.92_*  1.16    0.25   0.23_*
## 32  0.27  -0.44    -0.26    -0.80_*  1.04    0.22   0.16  
## 40 -0.12   0.20     0.12     0.36    1.22_*  0.03   0.16  
## 48 -0.35  -0.49     0.36    -0.55    0.80_*  0.22   0.05

Impacto dos Pontos Influentes

Observação 2

modelr4_sem2 = glm(trat ~ leuc + tmax, subset = -2, family = binomial)
impacto(modelr4, modelr4_sem2)[3:5]

##                impacto  pval_com_i  pval_sem_i
## (Intercept) -31.779726 0.021058970 0.010676263
## leuc         -1.303485 0.001251647 0.001803289
## tmax        -30.540927 0.017582847 0.009263493

A 5% de significância, removendo-se a observação 2, não há mudança inferencial no modelo.

Observação 32

modelr4_sem32 = glm(trat ~ leuc + tmax, subset = -32, family = binomial)
impacto(modelr4, modelr4_sem32)[3:5]

##               impacto  pval_com_i   pval_sem_i
## (Intercept)  11.05764 0.021058970 0.0449955910
## leuc        -15.84809 0.001251647 0.0009397309
## tmax         10.00241 0.017582847 0.0367816803

A 5% de significância, removendo-se a observação 32, não há mudança inferencial no modelo.

Observação 40

modelr4_sem40 = glm(trat ~ leuc + tmax, subset = -40, family = binomial)
impacto(modelr4, modelr4_sem40)[3:5]

##               impacto  pval_com_i  pval_sem_i
## (Intercept) -4.298702 0.021058970 0.017876120
## leuc         4.587500 0.001251647 0.002136859
## tmax        -3.948505 0.017582847 0.014971725

A 5% de significância, removendo-se a observação 40, não há mudança inferencial no modelo.

Observação 48

modelr4_sem48 = glm(trat ~ leuc + tmax, subset = -48, family = binomial)
impacto(modelr4, modelr4_sem48)[3:5]

##               impacto  pval_com_i  pval_sem_i
## (Intercept) -28.40658 0.021058970 0.012766838
## leuc        -28.35986 0.001251647 0.001124342
## tmax        -28.19501 0.017582847 0.010869361

A 5% de significância, removendo-se a observação 48, não há mudança inferencial no modelo.

2.6 Interpretação

Com essa análise, chegamos ao seguinte modelo para explicar a probabilidade de eficiência do tratamento:

\[ \log \left\{ \frac{\pi(\mathbf{x})}{1-\pi(\mathbf{x})} \right\} = 70.82 + 0.34 \cdot leuc - 0.07 \cdot tmax \; . \]

Como ilustração, ao fixarmos a temperatura máxima e considerarmos um aumento de 1% no número de células leucêmicas, a razão de chances de tratamento satisfatório, denotada por \(\psi_{CL}\), é estimada em

\[ \psi_I = e^{0.34} = 1.40 \; . \]

Isto significa que, sob uma mesma temperatura máxima, o aumento de 1% no número de células leucêmicas representa uma aumento de 40% na chance do tratamento ser satisfatório.

3. Taxa Anual de Câncer Nasal

3.1 Dados

Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 6, página 345). Dados disponíveis em: https://www.ime.usp.br/~giapaula/textoregressao.htm, arquivo canc1.dat.

Os dados a seguir são provenientes de um estudo de seguimento para estudar a associação entre a taxa anual de câncer nasal em trabalhadores de uma refinaria de níquel do País de Gales e algumas variáveis explicativas: idade no primeiro emprego (4 níveis), ano do primeiro emprego (4 níveis) e tempo decorrido desde o primeiro emprego (5 níveis). São também apresentados o número de casos de câncer nasal e o total de pessoas-anos para cada combinação desses três fatores. Nesse trabalho foi proposto um modelo log-linear com resposta de Poisson, sendo o número de casos de câncer nasal com offset dado por log(pessoas-anos).

As variáveis e suas classes:

• idade: 1 (<20), 2 (20-27), 3 (27.5 - 34.9), 4 (35+), em anos;

• ano: 1 (<1910), 2 (1910-1914), 3 (1915-1919), 4 (1920-1924), em anos;

• tempo: 1 (0-19), 2 (20-29), 3 (30- 39), 4 (40-49), 5 (50+), em anos.

Box-plots

3.2 Modelo Log-linear de Resposta Poisson

Denotamos por \(Y_{ijk}\) o número de casos de câncer nasal para o i-ésimo nível de idade no primeiro emprego, para o j-ésimo ano do primeiro emprego e para o k-ésimo tempo decorrido desde o primeiro emprego, \(i,j=1,...,4\) e \(k=1,...,5\). Supomos, então, que

\[ Y_{ijk} \sim Poisson(\lambda_{ijk} t_{ijk}), \]

em que \(\lambda_{ijk}\) é a taxa média de casos de câncer nasal por unidade de tempo para a idade i, ano j e tempo k. O modelo saturado, nesse caso, é dado por

\[ \log(\lambda_{ijk}t_{ijk}) = \alpha + \beta_i + \gamma_j + \delta_k, \]

em que \(\beta_i\) é o efeito da i-ésima classe de idade no primeiro emprego, sendo \(\beta_1 = 1(<20)\); \(\gamma_j\) o efeito da j-ésima classe de ano do primeiro emprego, sendo \(\gamma_1 = 1(<1910)\); \(\delta_k\) o efeito da k-ésima classe de tempo decorrido desde o primeiro emprego, sendo \(\delta_1 = 1(0-19)\).

Modelo Saturado

## 
## Call:
## glm(formula = ncasos ~ idade + ano + tempo + offset(log(total)), 
##     family = poisson)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7941  -0.7074  -0.3282   0.2736   2.4937  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -9.27178    1.31892  -7.030 2.07e-12 ***
## idade2       1.67276    0.75214   2.224  0.02615 *  
## idade3       2.48173    0.75909   3.269  0.00108 ** 
## idade4       3.42817    0.78160   4.386 1.15e-05 ***
## ano2         0.61895    0.37130   1.667  0.09552 .  
## ano3         0.05414    0.46808   0.116  0.90791    
## ano4        -1.12614    0.45294  -2.486  0.01291 *  
## tempo2       1.59815    1.04748   1.526  0.12708    
## tempo3       1.75121    1.05521   1.660  0.09700 .  
## tempo4       2.35486    1.07010   2.201  0.02776 *  
## tempo5       2.81758    1.11828   2.520  0.01175 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 135.681  on 71  degrees of freedom
## Residual deviance:  58.166  on 61  degrees of freedom
## AIC: 153.39
## 
## Number of Fisher Scoring iterations: 6

Envelope

3.3 Diagnóstico

Pelos gráficos acima, percebemos que as observações 62 e 63 se destacam nos gráficos de alavancagem e distância de Cook, caracterizando-se como pontos de alavanca e possíveis pontos de influência. A observação 52 se afasta um pouco do limite superior do intervalo de confiança para os resíduos componentes do desvio e será considerada como um ponto aberrante. A seguir, estudaremos a influência desses três pontos sobre o ajuste do modelo.

Pontos Influentes

Observação 52:

##                 impacto           pval_com_i           pval_sem_i
## (Intercept)   1.5128475 0.000000000002067477 0.000000000004356923
## idade2        3.6554899 0.026148616816039749 0.032007215605215911
## idade3        5.3753065 0.001078024818957924 0.001991338461317067
## idade4        2.4766547 0.000011541390822584 0.000018039837865180
## ano2          3.4807475 0.095523831910207815 0.108029138724031110
## ano3        272.8270147 0.907913016477262413 0.846449314051831769
## ano4         -1.8553932 0.012908715898225178 0.011467766786010202
## tempo2        0.2290205 0.127083241095111599 0.128300569536079490
## tempo3        0.6315893 0.096997436979855758 0.099518186750867349
## tempo4        0.9039952 0.027764429970676213 0.029392120989161854
## tempo5        7.7644850 0.011749976792630948 0.021946993370582788

Observação 62

##                 impacto           pval_com_i           pval_sem_i
## (Intercept)  -0.1251849 0.000000000002067477 0.000000000001868126
## idade2        2.1123292 0.026148616816039749 0.029302768207374805
## idade3        1.7230709 0.001078024818957924 0.001300403343998420
## idade4       -6.9946550 0.000011541390822584 0.000003030112064364
## ano2        -14.8348937 0.095523831910207815 0.054868569500244241
## ano3        218.9254636 0.907913016477262413 0.892509547081333787
## ano4         -7.5608046 0.012908715898225178 0.008225427383831935
## tempo2       -8.8173583 0.127083241095111599 0.096874999883543325
## tempo3        2.3170463 0.096997436979855758 0.105239733814455366
## tempo4       -0.2317686 0.027764429970676213 0.027524866857942212
## tempo5       -1.1706799 0.011749976792630948 0.010909475879866461

Observação 63:

##                 impacto           pval_com_i           pval_sem_i
## (Intercept)   1.5520620 0.000000000002067477 0.000000000004898349
## idade2       -1.5614680 0.026148616816039749 0.024150242028928428
## idade3       -0.6671066 0.001078024818957924 0.001021541390270312
## idade4       11.1067183 0.000011541390822584 0.000167202020944200
## ano2         28.7232661 0.095523831910207815 0.251882602591362581
## ano3        -67.1709199 0.907913016477262413 0.847118892595951078
## ano4         -0.4606752 0.012908715898225178 0.012672163470177361
## tempo2       -4.4846202 0.127083241095111599 0.110680742091644552
## tempo3       16.0756078 0.096997436979855758 0.168525091205826927
## tempo4        2.8124054 0.027764429970676213 0.032349724301361418
## tempo5        3.5454736 0.011749976792630948 0.014984580367553775

Ao nível de 5% de significância, nenhuma das remoções testadas acima ocasionou mudanças inferenciais nas estimativas dos parâmetros.

3.4 Interpretação

Com essa análise, confirmamos o modelo saturado para explicar a taxa anual de casos de câncer nasal.

\[ \log(\lambda_{ijk}t_{ijk}) = \alpha + \beta_i + \gamma_j + \delta_k \]

resumo = summary(model)
table = data.frame(Efeito = c("Constante","20-27","27.5-34.9","35+",
                              "1910-1914","1915-1919","1920-1924",
                              "20-29","30-39","40-49","50+"),
                   Parametros = c("alfa",
                                  "beta2","beta3","beta4",
                                  "gama2","gama3","gama4",
                                  "delta2","delta3","delta4","delta5"),
                   Estimativa = resumo$coefficients[,1],
                   E.Padrao = resumo$coefficients[,2])
knitr::kable(table, "pipe")

	Efeito	Parametros	Estimativa	E.Padrao
(Intercept)	Constante	alfa	-9.2717833	1.3189152
idade2	20-27	beta2	1.6727561	0.7521394
idade3	27.5-34.9	beta3	2.4817313	0.7590948
idade4	35+	beta4	3.4281689	0.7816028
ano2	1910-1914	gama2	0.6189461	0.3713040
ano3	1915-1919	gama3	0.0541431	0.4680769
ano4	1920-1924	gama4	-1.1261375	0.4529412
tempo2	20-29	delta2	1.5981513	1.0474835
tempo3	30-39	delta3	1.7512094	1.0552076
tempo4	40-49	delta4	2.3548632	1.0701010
tempo5	50+	delta5	2.8175829	1.1182812

Notamos, então, que as estimativas acima são significativamente diferentes de 0 e que há fortes indícios de que a taxa anual de câncer nasal cresce exponencialmente com o aumento da idade no primeiro emprego e o tempo decorrido desde o primeiro emprego, mas decresce exponencialmente com o aumento de anos do primeiro emprego.